FROM BROWNIAN DYNAMICS TO MARKOV CHAIN: 
AN ION CHANNEL EXAMPLE 
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Abstract. A discrete rate theory for general multi-ion channels is presented, in which the continuous dynamics 
of ion diffusion is reduced to transitions between Markovian discrete states. In an open channel, the ion perme- 
ation process involves three types of events: an ion entering the channel, an ion escaping from the channel, or an 
ion hopping between different energy minima in the channel. The continuous dynamics leads to a hierarchy of 
Fokker-Planck equations, indexed by channel occupancy. From these the mean escape times and splitting probabil- 
ities (denoting from which side an ion has escaped) can be calculated. By equating these with the corresponding 
expressions from the Markov model the Markovian transition rates can be determined. The theory is illustrated with 
a two-ion one-well channel. The stationary probability of states is compared with that from both Brownian dynamics 
simulation and the hierarchical Fokker-Planck equations. The conductivity of the channel is also studied, and the 
optimal geometry maximizing ion flux is computed. 
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1. Introduction. The membrane of a eukaryotic cell is mainly composed of a lipid 
bilayer, which is impermeable to water solvated ions |2|. Ion channels are nanopores formed 
by transmembrane proteins; they allow ions to flow through and act as biological valves 
connecting the intracellular with extracellular domains. Ion channels are the main mechanism 
by which cells control the intracellular concentration of chemical species, as well as the 
potential gradient across the membrane. As such they play important roles in maintaining 
various functions of plant, animal and human cells. 

There are two main features which distinguish ion channels from other nanoscale porous 
media. Firstly, they may be selective, distinguishing between the charge and size of ions; for 
example, the potassium K + channel conducts potassium ions at a rate 10 4 times faster than it 
does sodium ions 1101 . Secondly, their conformations may change between open and closed 
states in response to an external stimulus such as a voltage gradient, ligand binding, or pH 
value. 

The molecular structure of many ion channels has been revealed by X-ray crystallogra- 
phy in recent decades, which provides insight into their features and function. For example, 
the potassium K + channel is composed of four identical subunits which create a cavity con- 
necting the cell interior to a selectivity filter at the outer end of the pore iflOl . The narrow 
selectivity filter is only 12 A long and about 3 A wide, which forces potassium ions with 
Pauling radius 1 .33 A to shed their hydrating waters to enter and pass in a single-file fashion. 
The oxygen atoms of four carbonyl groups form four rings around the selectivity filter, which 
generate local minima called binding sites in the overall energy landscape to coordinate the 
dehydrated ions. 

Mathematical models for ion channels include molecular dynamics (MD), Brownian dy- 
namics (BD) and continuum theory (Poisson-Nernst-Planck equations) in descending order 
of resolution |[8]|7][T5]. Molecular dynamics provides the most detailed description by mim- 
icking the motions and interactions of all atoms (from membrane proteins to free ions and 
even individual water molecules) at the molecular level (3j [T8l [14) . Since the relaxation of 
water molecules happens at the fastest timescale of 1 fs, the time step of an MD simulation 
has to be very small, and one needs to evolve a system of thousands of particles up to of the 
order of 0.1 ms to observe ion conduction. Such a simulation is obviously computationally 
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intensive, but much shorter simulations (of the order of 10 ps) can be used to obtain informa- 
tion about the local potential energy and the effective diffusion coefficient of ions, which can 
then be fed into BD simulations. 

Brownian dynamics ll22l [9] [16] |6) is a more coarse-grained simulation which replaces 
the solvent molecules (water) as a continuum, and represent their influence by a dialectric 
constant and stochastic forcing. The fluctuations of membrane proteins are ignored and the 
channel is approximated by a solid boundary. Because the dynamics of water and proteins 
are no longer included, a relatively long time step can be used, which greatly reduces the 
computational cost. In this paper, we focus on this level of resolution, and introduce a discrete 
rate theory that is based on observations from BD. 

The continuum model |2T] [TTJ calculates the potential energy by a mean-field ap- 
proximation of average ion positions, which yields a Poisson equation, and then formulates 
a Boltzmann equation (in equilibrium) or a Nernst-Planck equation (in non-equilibrium) for 
the ion concentration. These continuum partial differential equations (PDEs) can be solved 
efficiently; however the individual ion-ion interaction is missing in this mean-field assump- 
tion which then fails to predict some properties (e.g. saturation). Comparisons of BD and 
continuum theories in different channel configurations are presented in l9l [T6ll . 

Recently several hybrid models combining MD and the theory of stochastic processes 
have been proposed, which are able to include molecular details and access long time scales 
while keeping computational cost low. One idea is to apply the Eyring rate theory to the ion 
permeation process using the potential of mean force (PMF) calculated using MD. This is 
based on the assumption that channels have some binding sites, and ions pass through by a 
hopping mechanism: an ion fluctuates around a certain site before it obtains enough energy 
to overcome the energy barrier and hops into the adjacent vacant site. This ion hopping 
mechanism has been revealed by MD in channels with binding sites J3j|4][l4). In addition, 
the single file diffusion constraint imposed by the narrowness of the channel assures that ions 
cannot cross each other in the channel. Therefore, the continuous dynamics of ion diffusion 
can be represented by transitions between discrete Markovian states. 

The Eyring rate theory was originally designed for chemical reactions in the 1930s, with 
transition rates proportional to the exponential of the energy barrier and distance between 
binding sites ifTTl (as shown in [ 8 1 this overestimates the physical barrier in the ion-crossing 
process). A novel theory was proposed recently in [T] for a one-dimensional channel with 
sawtooth-like PMF, in which the transition rates are not approximated using the energy barrier 
but are obtained as the product of total escape rate from one binding site and the splitting 
probability determining the relative chance of landing in each neighbouring site. H) showed 
that an optimal size of binding site maximizes the ionic flux if the applied voltage exceeds a 
threshold. They assume the channel is occupied by at most one ion, whereby the resulting 
system forms a single Markov chain, and the rates can be solved explicitly. In the multi-ion 
channel considered here, ion-ion interactions as well as the higher dimension of the energy 
landscape mean that the complexity of the rate theory is greatly increased. 

In this paper, we present a general discrete rate theory for a multi-ion channel, and com- 
pare it with BD. The ion permeation process involves ion hopping, ion escaping and ion 
entering. For the purposes of this work we assume ion entry rates are known and focus on 
calculating the other rates in terms of the mean escape time and splitting probability. Be- 
cause of the complicated network between states the rates are more intricately related to 
these quantities than in the single ion case. Moreover, since analytical solutions for the mean 
escape time and splitting probability are not available, these must be determined by solving 
the corresponding PDEs numerically. The theory is illustrated by a two-ion channel with one 
binding site and two ion sources. We show that, as with the one -ion channel, there exists an 
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FIG. 2.1. A schematic structure of a channel, x = —L and x = L are the (artificial) left and right bound- 
aries connecting large reservoirs of electrolyte in and outside the cell respectively. and £ represent the overall 
intracellular and extracellular environments, respectively. 



optimal shape for the external potential that allows a maximal flux. 

The structure of this paper is as follows. In Section [2] we introduce a general theory for 
a multi-ion channel with a maximal capacity of N ions. We first present BD simulations and 
formulate an equivalent cascade of hierarchical Fokker-Planck equations for the probability 
distribution of ions. An illustrative example of a 2-ion channel is discussed and the probability 
distribution from the histogram of BD and the solution of the Fokker-Planck equation are 
compared. Next, a discrete rate theory framework is presented in Section|3]and the transition 
rates calculated. The 2-ion channel is revisited in this framework, and the result is compared 
with that from BD. In Section |4] we apply the theory to study the dependence of channel 
conduction on different parameters such as the diffusion coefficient, ion entry rate and depth 
of potential wells. In particular, we study the effect of the geometry of the external potential 
in Section [5] We conclude by discussing the advantages and limitations of this method and 
possible applications and extensions in Section|6] 

2. Brownian dynamics. In this section, we present the theoretical framework of BD 
simulation. Since we are interested in studying the ion permeation process, which occurs on 
a time scale of 10~ 7 s, and since conformational changes occur on a timescale of 10~ 3 s, we 
assume that channel is always open and does not change its conformation. 

Since the channel is very narrow and the ions pass through in single file fl4l . we will 
suppose that the motion is one-dimensional, that is, the centres of the ions will be constrained 
to lie along a line. The generalisation to a fully three dimensional channel is algebraically 
complicated but conceptually straightforward. Since ions cannot pass each other in one di- 
mension, we may neglect the finite size of the ions and model them as point particles with 
charge. 

We define the maximal capacity of a channel to be N, so that it can hold up to N ions at 
one time. We denote the number of binding sites in the channel by M. The parameters N and 
M vary among different channels; for example, a germicidal A channel has two binding sites 
(M = 2) and single-ion occupation dominates (so that N = 1, or perhaps N — 2to allow for a 
knock-on effect) IT] [19] . 

At the ends of the channel the pore opens out into the intracellular and extracellular 
space. A full model would include (probably continuum) models of these spaces, which 
would then be joined (preferably matched in terms of matched asymptotic expansions, but 
more likely patched) to the channel model. For our present purposes we need to introduce 
(artificial) interfaces (i.e. points) at the left and right ends of the channel such that an ion 
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passing through these interfaces is taken to have left the channel and passed into the external 
domains. Without loss of generality, we suppose that the left interface connecting the channel 
to the intracellular domain lies at x = —L, and the right interface connecting the channel to 



the extracellular domain lies at x — L, as shown in Fig. 2.1 Thus an absorbing boundary 
condition is imposed at x — —L and x = L. 

In BD simulation of an ion channel the contribution of water molecules to the motion of 
a solute ion can be approximated by random collisions and an average frictional force in the 
evolution equation of the solute ion ll22l[T6l . The motion of a system of k ions is given by the 
Langevin equation 

mtdvi = -yvjdt +f^(x u ...,x k )dt + yV2DdWi, i = l,...,k, (2.1) 
where x, (f ) and v, (?) are the location and velocity of the i lh ion respectively. There are three 



forces on the right hand side of ( 2. 1 1. The first term corresponds to the frictional force exerted 
on the ion by averaging the effect of water molecules; y is the frictional drag coefficient, 
which depends on the surrounding fluid environment. Here we assume it to be uniform so 
that y is constant. The third term is the stochastic force generated by the random collisions 
of water molecules; W; is a Wiener process and D = kgT /y is the diffusion coefficient, where 
ks is the Boltzmann constant and T is the temperature. The second term fNx\,. . . ,Xk) is 
the overall electric force on the i™ ion, including interactions with all other k — 1 ions in the 
channel, fixed charges in the protein, and external field across the membrane. It depends on 
the locations of all ions, and can be obtained (along with the diffusion coefficient) from MD 
simulation. 

Note that a typical value of the diffusion coefficient in aqueous solutions at room temper- 
ature is D ~ 10~ 3 mm 2 s _1 , so that the ratio mj/y~ 10~ 14 s _1 . Since we usually take a time 
step At > 10~ 12 s in the simulation, the system is in an overdamped limit [8]. We may thus 
approximate p.l| l by the overdamped Langevin equation 

dx i = r —jf(x 1 ,... i x k )dt + ^2DdWi, i = l,...,k. (2.2) 



The boundary conditions on ( |2.2[ ) may be described as 

1. When the number of ions in the channel k is less than its capacity N, new ions 
are generated at the left (respectively right) end at a rate Hk (respectively Gk). In 
principle Hk and Gk depend on the current locations of the k ions in the channel 
x\ , . . . as well as the intracellular and extracellular environments J? and S. 
Since we are in the overdamped limit we cannot simply place the incoming ions at 
the ends of the channel: under Brownian motion they would immediately cross the 
boundary and leave the channel again. Instead we place them at a position within the 
channel given by the positional distribution function h{x;x\ ,x^) (or respectively 
g(x;xi ,Xk))- Note that h and g also depend on the positions of the existing ions. 
This is necessary since, for example, an ion entering the channel from the left must 
lie to the left of x\, while an ion entering from the right must lie to the right of x^. 
Thus, at the very least, h depends on x\ while g depends on x k . 

The functions h and g should be chosen to make the join with the outer model as 
smooth as possible, as in lfT2l [131 . Here we simply assume that h and g, and the 
rates Hk and Gk, are given. 

2. When X((t) < —L or xi(t) > L the 2 th ion is removed from the channel. 

3. If Xj(t ) > Xi+i(t) for some i then x, and jtj+i are switched. This enforces the single- 
file nature of the channel by preventing an ion overtaking its neighbour. This condi- 
tion is unlikely to occur with ions in a channel due to the strong Coulomb repulsion, 
but may be necessary if we are interested in neutral molecules. 
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Fig. 2.2. Hierarchical Fokker-Planck Equations describe the conservation of ions in the channel. For k-ion 
occupancy, the transitions to and from {k— \)-ion and (k + \)-ion occupancy by ions entering and escaping are 
demonstrated, along with internal transitions between states. 



2.1. Hierarchical Fokker-Planck equations. We denote by P k (x\,. . . ,x k ,t) the proba- 
bility density function for the event that there are k ions in the channel at positions xi,...,x k 
at time t . Since the number of ions in the channel may run from zero to the channel capacity 
N, we have 1 such probability density functions. The probability of no ion in the channel 
(i.e. k = 0) is denoted by Po(t), and is independent of the spatial variable. We label the ions 
by the order of their locations, such that x, < xj for i < j. Then the stochastic process ( |2.2| > is 
equivalent to the following hierarchical system of Fokker-Planck equations: 



d t P k (x u ...,x k ,t)=DV ■ (vP k (xi , . . . ,x k ,t) - P k (xi , . . . ,x k ,t) j^^k(xi , . . . ,x k ,t) 

- \ H k (x\ , . . . ,Xk) + G k (x\ ,x k )jP k (xi , . . . ,x k ,t) 

+ Hk-i(*2, ■ ■ ■ ,x k )h(x l ;x 2 , ■ ■ ■ ,x k )P k _i(x 2 , . . . ,x k ,t) 
+ G k ^i(xi ) ...,x k -i)g(x k ;x 1 ,...,x k - 1 )P k -i(xi,...,x k -i,t) 
+ T k+ i (xi,... ,x k )+R k+ \ [x\,... ,x k ), 



(2.3a) 



where V = (d X} , . . .,d Xk ), F, - (/* . . .,/*) e R k , and 



T k +i(x2, ■ ■ ■ ,x M ) = D {-j^~ - p k+i ~j~ff\ + X ^j (- L > x 2> ■ ■ ■ >**+l)) 
R k+ l(x\,...,x k ) = -D^^±i -P k+i ^f^l^j(x u ...,x k ,L), 

where k — 0, . . . ,N and we use the convention that P_i = Pn+i = 0. Since only one ion can 
escape or enter at any one time, P k is coupled only to the neighbouring states P k _\ and P k +\. 
Note that Hpj = Gn = 0, since no ions can enter when the channel is fully occupied. 



The first two terms (i.e. the first line) on the right-hand side of (2.3a i correspond to ion 



diffusion and ion drift respectively, where the drift term includes the external potential as well 
as ion-ion interactions. The third term corresponds to a new ion entering the A:-ion channel 
from intracellular or extracellular solution; this term is negative since such an event leads to a 
transition from a A:-ion channel to a (k+ l)-ion channel. The fourth and fifth terms correspond 
to a new ion entering a (k — l)-ion channel from the left and right respectively. The sixth and 
seventh terms (i.e. the last line) of (2.3ai correspond to ions leaving a (k+ l)-ion channel 
from the left and right respectively. 
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The boundary conditions on (2.3ai are 



P k (-L,x 2 ,...,x k )=0, (2.3b) 
P k (xu...,x k -i,L)=0, (2.3c) 



along with the no-flux condition on the interface x,- = x, 



i+U 



lim (^- p ^) = lim (p—^rr^ ,13d) 



for i = 1, . . . , k — 1, which ensures that the ions are correctly labelled. The reason we have 
had to write this as a limit is that the inter-ion potential tends to infinity as xi — > Jt,-+i, while P k 
tends to zero. A local analysis shows that we need P k to tend to zero faster than (xj+i — x;) 2 . 
We note the following normalisation condition, which holds at all times t , 

N , 

Pk(xi,...,Xk,t)dxi---dxk = l, (2.4) 

*=<r r * 

where is the available state space when there are k ions in the channel, namely = 
{(x\,. . . ,Xk) ■ x\ < *2 < • • • < x k}- We will usually be interested in the steady state; in that 
case we solve the coupled hierarchical Fokker-Planck equations for the stationary probability 
distribution = lim r ^< x ,Pj : (f ) for k = 0, 1, . . . ,N. 

2.2. An example with N = 2. We exemplify the theory above with a simple channel 
that is selective to cations with elementary charge e = 1.6 x 10~ 19 C. The selectivity of this 
type of channel is generally caused by negative charged boundary proteins, which decrease 
the energy barrier imposed by the narrow structure and assist the permeation of cations. For 
example, the oxygen atoms of four carbonyl groups in the selectivity filter of the potassium 
channel can be modelled by putting four negative partial charges equally spaced on a ring of 
radius d that is perpendicular to the x-axis J9)- 

We consider the simplest possible example of multi-ion channel with capacity N = 2 
and a single binding site M = 1. The binding site is located at the position x = % and is 
a potential well generated by a ring of fixed partial negative charges a distance d from the 
channel axis. By Coulomb's law, the potential energy <&i(jci) seen by one cation at x\ with 
charge e traversing through the channel is 

*i(*i) = A( Tl ~ k t o +^O- ( 2 - 5a ) 

where k e the Coulomb force constant and Z the total fixed charge on the ring, and U is the 
constant field, which imposes a potential difference 2UL across the channel [— L,L]. This 
potential difference is small compared to the potential well, and does not change the shape of 
the potential well but merely tilts it by a small angle. The force on the ion due to the potential 

is 

fl = -k B T-±. 

ax\ 

When there are two cations in the channel, at positions x\ and x%, the overall potential 
energy < t>2(xi,x : 2), including the interaction between the two free ions, is 

*2(*i,xa) = = J7 + = J9 + +U(x l+ x 2 )). (2.5b) 

k B T \^(xi-%) 2 +d 2 v / (x 2 -^) 2 + d 2 \xi-xz\ J 
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The forces on the two ions are then 



Finally we need to specify the entry rates and G^ and entry distribution functions h and 
g. We choose the simplest possible model for the entry distribution function. We suppose 
that the ions entering from the left are all placed at a position x_ near the left-hand end of the 
channel, while ions enterying from the right are placed at a position x + near the right-hand 
end of the channel, that is 

h{x) = 8(x — X-), g(x) = 8(x — x+). 

We have to be careful in implementing this condition that we preserve the order of the ions in 
the channel. We choose to do this as follows: if we are attempting to place an ion at position 
and the position of the existing ion x\ < x~, then we abandon the insertion of the new ion. 
A similar procedure is implemented at the right-hand end. In effect this means that the rate 
of entry is chosen to be zero whenever the position x\ of the existing ion is such that x\ < x_ 
or x\ > x + . (An alternative procedure would be to modify the distribution functions h and g 
so that h = if X\ < x_ and g = if X\ > x + , but this would mean altering them from the 
present 5-functions.) 

In general the entry rates may be functions of the current ion numbers and locations as 
well as the intracellular J? and extracellular $ environments. However, for this illustrative 
example we suppose that they are constant subject to the constraint set out above. Thus we 
choose 

Hq = X, G(,=pL, H\ = X&{x\ — jc_), Gi = pi®{x + — xi), 

where © is the Heaviside function. Recall that H2 — G2 — O since the channel is then fully 
occupied. To run the Brownian simulation, we set the time step At = 100 ns, and the physical 
parameters as 

L = lnm, jc± = ±0.9nm, ^ = Onm, c/ = 0.5nm, D = lnm 2 • ns~' , A=/i=5ns~ 1 , 
T=298K, k B = 1.38 x lO^J-KT 1 , t/ = 0V-nrrr\ Z = e. (2.6) 

We use the nanometer as the unit of length and the nanosecond as the unit of time. We evolve 
( |2.2| > for 2 x 10 9 iterations until a dynamic equilibrium is reached. During the simulation the 
number of ions in the channel varies in time as ions enter and leave. We record the number 
of ions and their locations at each time step. We find that the proportion of time spent with k 
ions in the channel, Jk say, is given by 

7 « 0.0000, 7i« 0.8986, 7 2 ~ 0.1014. (2.7) 

Thus, for these parameters, the channel is almost never empty, and for nearly 90% of the 
time there is just one ion in the channel, with two ions the remaining 10% of the time. The 



histograms of 2-ion distribution and 1-ion distribution are plotted in Fig. 2.3(a) and Fig. 



2.4(a) respectively. 

The stationary probability distributions ^2(^1, ^2), P\{x\) and Pq satisfy the stationary 
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FIG. 2.3. Using the parameters in (2j6j, the stationary probability density of a 2-ion channel is computed as 
(a) histogram from Brownian dynamics simulation; (b) solution of j2.8a) - (2~8e) . Here x\ is the position of the first 
ion andx2 is the position of the second ion; since we label the ions such thatx\ < X2 the state space is a triangle. 



Fokker-Planck equations for a two-ion channel, namely 

= DV-(vP 2 (x 1 ,x 2 )+P 2 {x l ,x 2 )V<t> 2 {x h x 2 )' S } 

+ X&(x 2 —X-) 8{x\ -x_)Py{x 2 ) +jU©(x + — x\)8(x2— x + )Pi(x\), (2.8a) 

0=D^l^(xi)+P\(xi)^(x 1 ))-(?i&(x 1 -x-)+^®(x + -xi))Pi(xi) 



0: 



h A 8 {xi -x-)P + n8 (x\ - x+)Po 

( dP 2 ~ d<P 2 \ , ( dP 2 ~ d<Z> 2 \ , 

ID j± +Pi^ (-L, Xl )-D _* (x u L), 
\ ax\ dx\ J \ dx 2 dx 2 I 

( dPi ~d4>i\, , ( dP x ~d4>i\ /N 



with the boundary conditions 

P 2 {-L,x 2 )=P 2 {x u L)=0 } 

and 



Pi (-L) =P 1 (L) =0. 



lim [f*+ft^ 



= lim 



dP> + ^ 5*2 

<9x2 dx 2 



(2.8b) 
(2.8c) 

(2.8d) 
(2.8e) 



We solve ( |2.8a[ )-( |2~8e| ) by the finite element PDE solver C ortisol with 28800 elements. 
The stationary distribution P 2 {x\,x 2 ) is shown in Fig. 2.3(b) and P\(x\) is shown in Fig. 



2.4(b) We see that these agree with the histograms in Fig. 2.3(a) and Fig. 2.4(a) obtained 



from Brownian dynamics simulations. 

We see that P2(xi,X2) is localised around two discrete states near (x_,<!;) and (^,x + ), 
while P\(x\) is localised around x\ = ^. The most likely path between the two states of 
P 2 (xi,x 2 ) can also be faintly seen. 

This localisation of P 2 and Pi motivates the definition of a small number of discrete states 
which the system can adopt, which is the basis for the discrete transition rate theory described 
in the next section. 




3. Discrete transition rate theory. We saw in our two-ion example (Figs. 2.3 and 



2.4 1 that P2 was mainly localised around two regions in state space, while Pi was mainly 
localised around one region. Suppose, in general, that when there are k ions in the channel the 
stationary probability distribution Pk{x\ ,Xk) is mainly localised around small regions. 
Let us denote these regions by 

Sfcr, for* = l,...,Z*; 
then Pk is very small outside U^S^\ so that 

The idea of discrete rate theory is to replace the continuous variable P^ with a set of discrete 
probabilities corresponding to the states so that 

H = / (i) Pk(x\,...,x k )dxi--- dx k 

Js k 

is the (stationary) probability that (*!>••• ,*k) & Note that pj^ is just a number: it is 
independent of spatial variables. In total there are L% = Y!H=o^k states in the channel, and the 
sum of the probabilities of all L-^ states is unity according to ( |2.4) , that is, 

k=0i=l 

We now imagine a Markov chain in which the channel undergoes transitions from one of 
these discrete states to another, with the transition probabilities dependent only on the current 
state (i.e. no past history is involved). This Markov chain is illustrated in Fig. |2.2| Such 
Markov chains for ion channels have been previously considered for a single ion in a many- 
well channel H). However, multiple occupancy of the channel leads to a more complicated 
transition structure. 
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Since only one ion can enter or leave at once (so that P k is coupled only to P k -\ and Pk+i) 
we see that may have transitions to and from only the states sjp , and s[jj_ r The 
general master equation for the time-dependent probability pjp (t) is of the form 




influx 



E/r , +Er ,, Kw, (3.i) 



Outflux 

where Gf^ ^ is the transition rate from to S^L, JS^ ^ is the transition rate from to 

Su_ v and 7^' J ' is the transition rate from to Sjp. Thus a describes the influx of a new 
ion, j3 describes the loss of an ion to the intracellular or extracellular environment, and y 
describes a hopping of the ions within the channel. In fact we expect many of these rates to 
be zero, since, for example, when we add a new ion to a channel it must occupy either the 
left-most or rightmost potential well. 

The entry rates for new ions (xj,'^ may be determined from H k and G k , which for the 
present purposes we are assuming are given. The ion escape rates j3t and hopping rates 

i k ,j ' can be computed from the notation of mean escape time and splitting probability, as 
described below. 

To define the mean escape time we set all the influx probabilities to zero. We then 
suppose that the channel initially contains k ions located at positions x\, . . . ,x k . We define 
the mean escape time %^{x\,... ,x k ) to be the average time before the channel undergoes a 
transition to a (k— l)-ion configuration, that is, the average time for one ion to leave the 
channel. Using the backward-Kolmogorov equation [20] it can be shown that satisfies 

At* - V<D r Vt* = ~, (x u ...,x k )e r k (3.2a) 
% k = if x\ = —L or x k = L. (3.2b) 
Then the mean escape time from state S$ is given by 

[ { , ) P k dx 1 ---dx k 

3 k 

We now determine a similar expression for T^fsP] using the discrete transition rate model. 
Equating the two expressions will then provide information on the rates JS^ and % ■ 

To this end suppose that the channel is initially in the state si' , so that P k \o) = 1, 
Pm 1 (0) = otherwise. As before the influx rates a k are set to be zero. The master equation 



this solution we can determine the mean escape time Tfc[SjP] as 



( |3 . 1 ] > for P k then decouples from those for P k \ and P^K, and we can solve for P k . Given 

f 

^[Sf] = '_ } _ f ~ — • (3.4) 



EEf AWw* 
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In calculating the mean escape time we have not distinguished between the case that the first 
ion leaves from left end into intracellular electrolyte ^ and the case where the last ion leaves 
from right end into extracellular electrolyte $ . However, it is important that the discrete state 
model gets the ratio of these probabilities correct, since this is what causes a net ionic flux 
through the channel. Thus the second piece of information we use to determine the rates 
j3j and yi is the splitting probability Pk(xi,. . . This is defined to be the probability 
of the first ion to exit was x\ from the left-hand side of the channel, under the condition that 
an ion-escaping event from a Ai-ion to a (£— l)-ion channel has occurred, given that the k ions 
started in positions (x\ , . . . initially. The splitting probability function p k satisfies, 

Ap k -V<P k -Vp k = for (*!,...,**) er* (3.5a) 
p k = 1 on x\ = —L. pit = on x k = L. (3.5b) 

As with T k , we can now calculate the splitting probability for state S$ as 

j{i)PkPkAx\---Axk 

Pk[sf] = a—- ^ 

J j,)P k dxi ■■■6x k 

3 k 

To calculate the splitting probability from the Markov chain we need to separate /V into 
two individual rates representing the case that an ion leaves to the right into the extracellular 
domain, and the case than an ion moves to the left into the intracellular domain, that is, we 
write 

Then the probability that an ion escapes to the left given that it escapes is 

Pk [s£°] = ?= — ?= ■ (3-7) 

LL () Pk iLj) p { k j) m+LL L PpMpPw 

I j JO i j JO 

Note that, as in the case of the escape time Tk, the right-hand side depends on sf through the 
initial condition on . 



By equating ( 3.3 i with ( |3.4] > and ( 3.6 1 with ( 3.7 1 we have a number of equations to help 
determine the unknown rates pj: and "K . Since only the left-most (respectively right- 
most) ion can escape from the left-hand side of the channel (respectively right-hand side), 
many of the rates j3j[ will infact be zero. If we still do not have enough equations to 

determine the remaining j3j and jk'^ , then it will be necessary to determine some of the 
transition rates between internal states. Since these do not involve a change in the number of 
ions in the channel, they may be determined by standard techniques. 

Note that to determine the net flow of ions through the channel we will also have to 
distingiush between ion entry from the left and from the right, that is, we should also split 

a p = a -( i J) + a+ (iV) . 

However, in most cases (at least) one of these rates will be zero, since it is not possible to 
have the same transition between two states occuring with an ion entering from either side. 
The one case where this is possible is the transition between an empty channel and a one-ion 
channel, which occurs in our example below. 
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FIG. 3.1. The transitions between the four different states: the three circles represent the left entry point, the 
binding site and the right entry point; a (green) filled circle indicates the presence of an ion. 



3.1. Example of two-ion channel. Now we revisit the example of two-ion channel in 
Section [272] and illustrate the rate theory using the parameters in \2.6\ . 

We have seen that the channel can exist in a 2-ion, 1-ion or 0-ion state. From Fig. |2.4| 
we see that P\(x{) is localised around the single regi on x \ = so that there is only one 



metastable state with one ion in the channel. From Fig. 2.3 we see that P2(xi ,^2) is localised 
around the two states and (^,x + ). Thus there are two metastable states with two ions 

in the channel. Thus our Markov chain comprises the four states 

:{(*-,$)}. 4 2) :{(<^+)}, S? ] :{Sh S^:{}. (3.8) 

Thus L2 = 2, L\ = 1, Lq = 1 and overall there are Lz — 4 states for this channel. These states, 
and the transitions between them, are illustrated in Fig. |3.1| The circle at center represents 
the binding site x = £, , and two other circles represent the left and right entry positions x=x±. 
A (green) filled circle represents a position occupied by an ion. Note that for the transitions 
between si and Sj it is important to distinguish between ions entering and leaving from 
the right and from the left, so that we can calculate the net flow of ions through the channel. 
Let us first consider the ion entry rates. We find 

o+ (M) = 0, ar (M) =A, a+^/i. <h™ = 0, a +(11) = M , a Q (L1) =X. 

Note that the two zero values arise because the transition from to occurs via an ion 
entering from the left, while that from Sj to occurs via an ion entering from the right. 
Note also that Ot^ = f° r a ^ h j smce with two ions in the channel is already full to capacity. 

Let us now consider the ion leaving rates j3t . In principle we have six of these to 
determine. However, since the transition from to 5j must occur via an ion leaving from 
the left we know that jS^ = 0. Similarly the transition from 5^ to must occur via an 
ion leaving from the right, so we know that /3 9 ^' 2 ^ = 0. This leaves /3 2 , > Pi 

and j3j 1 to determine. To these we must add the two hopping rates J^' 2 ^ an< ^ • 

Denoting the state of the system by the probability vector P = {P^\p^\p^ l ',p{p) T , the 
master equation governing the evolution of the Markov chain is then 

dP 

— = ^-P(t), (3.9a) 
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where the 4x4 transition matrix & is given by 



(1,2) 



r 2 



.(2,1) 



72 



-/3 2 +(1 ' 2) - 



-y 2 



1.2) 



-(i,i) 



h(2,l) 



V 



-(1,1) 







?+(l,2) 







-(1,1) 



K2,l) 



_£+(i,i)_£-(M) 
^(i,D +j8r (i.i) 



-a +{l ' l) -a- {1 ' l) 



(3.9b) 



As expected, the sum of each column of the matrix 3F is zero (since the system ( 3.9a i con- 
serves probability), so the matrix is rank deficient. The stationary probability P is the eigen- 
vector associated with zero eigenvalue of matrix ST , 

To calculate the mean escape time and splitting probability we set the all entry rates 
to zero and solve ( 3.9a ». To emphasize that this is an auxilliary problem and not the true 

Markov chain we denote the probability of lying in each state by ?2 (*)i?2 (0> ^i^Mj^o^ W 
respectively. Then ( 3.9a| > is 



df 

A (2) 

df 



r 2 



= - ft 



.(2,1) 



"(1,2) 



(1,2) \ (2) . ,,(2,1) H) 



^2+72 ' ^2 



(i) 



df 



(i) 



df 



!ft -(M) s (t) +ft+ a V . 



,fM +ft -W)),f). 



(1,1) 



The first two equations decouple. We start by considering the state that is, we solve 
(3.10ai-(3.10b i subject to the initial conditions #2 (0) = 1 an d <7 2 2 '(0) = 0- This gives 



(3.10a) 
(3.10b) 
(3.10c) 
(3.10d) 



4 2) 



1 



+ (1,2) , ,,(1,2) 



X\ — A 2 

1 



J2,l) 



7^ 



exp(Ai f ) 



A 2 +j8 2 +(1 ' 2) +ri 1 ' 2) 



A2 — A] 1 

where Ai , A 2 are two eigenvalues satisfying 

A 1+ A2 = -fft 2 -^+ft 2 + ^- 



r 2 



.(2,1) 



exp(A 2 0) 



(3.11) 



,,(2,1 



(1,2) 



A 1 A2=ft2 J(1 ' 1) ft 2 +(1,2) +^ (1 ' 1 VP ) +j3 2 +(1 ' 2) 7f 1) . 
(t)i 



Using (3.4 1, the mean escape time T 2 [S 2 ] is 



ftT^^W+ftT^W* 
fe +(1 ' 2) + yf 1) + y^ 2) 
^ +(1 - 2) 7 f 1) +^ (U) 7^ 2) +^ (U) /3 2 +(1 ' 2) 
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(3.12a) 



and, using (3.7 I, the splitting probability p2\$2 ] 1& 



p 2 [4 !) 



q\ \t)At 



ft- (I,I) ^ 1} (0+fe +(I ^W* 



-(1.1) (1,2) 
12 



-(1,1) (1,2) 
/2 



Similarly, by applying the initial conditions q^p(0) = and q^ifi) = 1, we find 



i-p 2 [4 2) ] = 



r 2 



1.2) 



,+(1,2) ,,(2,1) 



(l,l)o+(l,2) 



J8- (M) J8 2 +(1 ' 2) - 



^ (u) j3 2 +(1 ' 2) +i3- (u Vl 1 ' 2) +i3 2 +(1 ' 2) 7f' 1) ' 



(3.12b) 



(3.12c) 
(3.12d) 



Finally we have to consider the escape time and splitting probability for state S\ . With the 



initial condition q 2 l> (0) = q^ify — 9n^(^) = ^, #1^(0) = 1' equation ( 3. 10c I decouples and 
is easily solved to give 



,(1) 



tJS 



tqf\t)dt 
q[ l) dt l 

, r (l.l) +/s +(l.l)- 



,-(1,1) 



(1,1)' 



(3.12e) 



(3.12f) 



The mean escape times T2, T[ and the splitting probability P2, pi can be obtained by solving 
p. 2} and (p3J, respectively with k = 2, 1 by Comsol. Fixing £/ = 0, <!; = 0, the mean escape 



time T2 in triangular domain T2 is plotted in Fig. 3.2(a) the splitting probability P2 in T2 
is plotted in Fig. 3.2(b) Since there is no applied field across the channel (U — 0), and the 
potential well is at the center, the external potential is symmetric with x = 0, therefore both 
functions are symmetric with x\ +.12 = 0. We take the values of functions at the centre of the 
different states and obtain 



T 2 [S ( 2 l) ] = T 2 (jc_ ,| ) = 0.01 113, p 2 [Si 1 '] = p 2 (x- , 4 ) = 0.9964, 



(1)1 



x 2 ^,x + ) =0.01113, p 2 [sf }=p 2 {$,x 



Ti [S^] = ti (4 ) = 3.9385 xlO 3 , P1 [S[ 



= 0.0036, 
Pl (4)=0.5. 



Equating these expressions to those of ( |3.12| i we find six equations for the six unknown rates. 
Solving these we find 



j8 2 _(M) = j8 2 +(1,2) = 89.8283, ff l > = f 2 ^> = 0.3361, p[ 



,(2,1) 



(1.2) 



1,1 



= P 



= 1.2695 x 10 



-4 



Note that by using the mean escape time and the splitting probability we have not had to 
estimate the internal hopping rates y 2 from the Fokker-Planck equation, but have been able 
to determine them from the auxiliary problems we have solved. We will see in Section[5]that 
this is especially useful when the internal states are not so well defined. 
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(a) T2(xi,X2) 



(b) P2{X\,X 2 ) 



FIG. 3.2. The parameters are given in (2.6) . (a) Mean escape time Vi(x\ ,x%) for 2-ion transiting into \-ion 
state, (b) Splitting probability f>z (jci , x% ) of ion exiting from left side under the condition that an ion escaping event 
occurs. 



Since there is no external potential gradient in this case (U — 0) the rates are symmetric, 
and there is no net flux through the channel. However, we can already make some observa- 
tions. Firstly, the rates jSj^ are tiny compared to the others. Thus the channel will switch 
between single and double occupancy, but will almost never be empty of ions. We will con- 
firm this when we consider the equilibrium occupancy of the channel in the next section. 



Secondly, the exit rates j3 2 and jS^^'^ are about 270 times as large as the hopping rates 

(2 1) (1 2) 

72 ' and y-J ' and. This means that, for these values of the parameters, an incoming ion 
enters and leaves about 270 times before it manages to replace the bound ion in the potential 
well at the centre of the channel. 

3.2. Stationary probability of each state. Now that we have determined all the transi- 
tion rates in Fig. |3.1| the stationary probability for the number of ions in the channel can be 
calculated explicitly as 



(1,2) 



J(2) 



J(l) 



a, 



^2[4 



-i,i 



KM) 



a 



1 + - 



i 



CM) ,„-KM) 



a 



-(i,i) 



a, 



a, 



-(2,1) 



-(2,1) 



(2), 



(2), 



: 0.1002, 



1 0.8998, 



J(i) 



Oh ^MM) 



KM) 



1 



-(i,i) +(i,i) 



a 



-(i,i) 



a 



+ (2,1) 



i 2 h 



; 2.285 x 10 



-5 



(3.13) 

We see that these agree well with the probabilities obtained from a Brownian dynamics simu- 



lation in ( 2.7 1. The same probabilities may be obtained by integrating the solutions of Fokker- 



Planck equations over the configuration space, which gives 



h = 
h = 



P2(xi,x 2 )dxidx 2 tt 0.1001, 

-LJ-L 
L _ 

Pi(xi)dxi 



1 0.8991, 7 = Po~ 8.3x10 



-4 



(3.14) 
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4. Current-voltage curve. The most important characteristic of an ion channel is its 
conductance. In this section, we investigate the channel conductance by examining the 
current-voltage curve for various values of the channel parameters. 

The potential difference across a channel is usually around 100 ~ 200 mV. We therefore 
vary the potential gradient U in the range [—0.1,0.1] V nm _1 , which gives a voltage drop in 
the range [—200,200] mV for a channel of length 2 nm. Following the framework in Section 
[3] we compute the transition rates corresponding to each given value of U . 



By examining the transition network in Fig. 3. 1 we see that there are two different paths 
which lead to ions moving from the intracellular (left-hand side) to the extracellular (right- 
hand side) domain, namely 

-(1,1) (2.1) d+(1.2) o-(l.l) -(1.1) 

PATH 1 : sf)^)i->^*-fSf) PATH 2 : S^^S^^sfK 

Both paths start with a channel with one ion bound at the potential well. In Path 1 another ion 
first enters the channel from the left-hand source to produce a two-ion channel. The two ions 
then hop to the right so the new ion lies in the potential well at the centre of the channel. The 
ion released from this well then exits the channel at the right. Thus in Path 1 we can think of 
a new ion coming in and knocking the present ion out the other side. In Path 2 the ion in the 
channel first leaves from the right to leave an empty channel, and then a new ion enters from 
the left. 

By considering the transition rates we can determine the relative importance of each of 



these mechanisms. For the parameters in ( 2.6 » the rate j8 t is tiny and Path 1 dominates 



the current. Note that at equilibrium the ion flux entering from the left is balanced by the flux 
which leaves from right for each path, so that 

„-(U)p(l) _ „+(l,2)p(2) „-(Ll)n(I) _ o-(U)p(l) 

a l 1 — h>2 r 2 ' u ~~ Pi 1 

Similarly ions flow from right to left via the paths, 

+(2,1) (1.2) d-(I-I) o+(l-l) +(1.1) 

PATH 3 : S« , PATH 4 : — > — > s[ X \ 

and in equilibrium 

+ (2,1)~(1) _ n-(l,l)p(l) „ + (l-l)p(l) _ o+(U)p(l) 

"1 1 — P2 r 2 > a — Pi 1 

Combining the current from each path the net current is given by 

•.(fr™^ + ti {L2) P 2 m - a^P^ - a+^/f ) , (4.1) 

Next we study how the conductance of the channel varies with the external potential ene rgy 

the 



■e 



4.1 



parameter d/L and the dimensionless entry rates L 2 X/D and fj,L 2 /D. We plot in Fig 
current-voltage (I-V) curve for different entry rates when d/L = 0.5. For large voltages the 
current saturates since it is limited by the entry rate of ions A and jj,; this effect is illustrated 



in Fig. 4.1(b) In Fig. 4.2 we plot the I-V curve for various values of d with a fixed entry 
rate. The slopes of these curves give the conductance of channel, which grows initially with 
increasing voltage until diminishing as saturation sets in. We see that as the dimensionless 
entry rate increases, or the potential well gets shallower, the conductance of the channel 
increases. 
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FIG. 4.1. I-VcurveforL 2 ?i/D = L 2 fl/D=l,5, 10 with d/L = 0.5; all other parameters are as in (2~6) . (a) 

for small voltages, for which internal transitions in the channel are rate limiting; (b) for larger voltages, at which 
saturation occurs due to the finite rate of entry of ions from the hulk. 
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FIG. 4.2. I-V curve for d/L = 0.3, 0.4, 0.5 withL 2 X/D = L 2 n/D = 5. All other parameters are as in (2.6) . 



5. Optimal geometry of potential. We showed that a channel with shallower potential 
well has larger conductivity in Fig. |4.2| if all other physical parameters are fixed. However, 
the depth of the potential well is not the only factor that determines the ion flux. Abad et al 
[l] studied the flux through a channel of capacity N = 1 with symmetric M-shape potential 
energy, and showed that there exists a critial ratio a of the width of potential well over the 
length of channel, at which the flux is maximized. This optimal geometry of potential energy 
requires the potential well is neither too narrow (a — > 0) nor too wide (a — > 1). 

We now perform a similar analysis of the multi-ion channel. For our N = 2 example we 
study how the shape of an external potential on the protein boundary affects the conductivity 
of the channel. We vary the distance d and carefully choose the external charge as 

2-(l+0.5 2 r°' 5 



d~ l L-(l +d 2 L~ 

so that the depths of the resulting potential wells are all the same (and so that Z = e, when 



d/L = 0.5); the variation of the potential well with d is shown in Fig. 5.1(b) Th e resul ting 
group of potential wells with applied field U = —0.05 V run -1 are plotted in Fig. 5.1(a) the 
darker curve corresponds to Z = e, d = 0.5L, the dotted curve shows how the potential wells 
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(a) <I> vs. x 



0.4 0.6 

d (nm) 

(b) Zvs. d 



FIG. 5.1. (a) A group of potential energy wells for d = 0.1,0.2,0.3, . . . ,0.9, 1.0 run are plotted by solid curves 
respectively from narrow to wide well, the darker curve corresponds to Z = e, d/L = 0.5, and the dotted curve 
corresponds the energy drop due to applied field U = —0.05 Vnm~K (b) The external charge Z as a function of d. 



tilt with the applied field. Obviously, when d is small, we have a steep potential drop near the 
binding site, and a relatively flat energy landscape near the ends of the channel. When d/L 
gets large, the potential well tends to a curve that is steeper near the end of the channel and 
flatter near the binding site. 

Now we fix U = -0.05 V nm" 1 , 



5.2 



the stationary 



k = /J = 5 ns -1 , and plot in Fig. 
probability density function Pi(x\,X2) obtained by solving stationary Fokker-Planck equa- 
tions numerically for various values of d. Recall that since x\ < x%, we only look at the upper 
left triangular domain. We observe that when the potential well is very narrow and steep at 
the binding site with d = 0.1 nm, the probability distribution of two ions is localized at two 
tiny spots around (jc_,£) and (£,,x + ). As d increases until d = 0.5 nm, the two spots grow a 
little larger and shift slightly away from (jc_,§) and (£,jt+). Thus for d £ [0.1,0.5] nm we 
can justify the use of our four-state rate theory. For d between 0.5 nm and 0.6 nm, a ridge of 
high probability distribution emerges that connects previous two spots at its two tails. This 
implies that instead of the two ions being trapped at one of two states, the two can wander 
back and forth freely between these two states. As d increases even further to d = 0.9 nm, the 
ridge shrinks to its center peak, which corresponds to the two ions both sitting in the potential 
well. Thus, for larger values of d, we have effectively only one state for the two-ion occupied 
channel. 

Thus, for large d, we can define a single-chained three-state system 



S< 1 ':{((*- + S)A($+*+)/2)}, 



: {}, 



(5.1) 



which is similar to four-state system in Fig. 



3.1 



except that si — ^2' are combined and 



(2) 



there is no hopping y^' 2 ^ and y^ 2 ' 1 ' between them. Following the framework in Section [i] 



all the entry rates a ( 



-(1,1) 



a 



-(1,1) 



jj. and a, 



+(1,1) 



a, 



escaping rates are easily calculated as 



!-(l,l) 



,+(1,1) 



-(1,1) 



X are prescribed, and the 



1,2. 



In Fig. |5.3(a)| we plot the mean escape time T2 from the left state 5^ = (—0.9, 0) (black solid 
(2) 

curve) and the right state S 2 = (0,0.9) (black dash-dotted curve) in four-state formulation, 
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(a) i = 0.1 nm 



(b) d = 0.2 nm 



(c) d = 0.3 nm 
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FIG. 5.2. Fixing X = fi = 5 bj -1 , {7 = -0.05 V nm , we pZof ?/ie stationary probability distribution Pi(x,y) 
for d = 0.1,0.2.0.3, 0.4,0.5,0.6, 0.7,0.8,0.9 nm, which are obtained by solving (2JSJ numerically. It shows the 
change from two distinct states at (—0.9,0) and (0,0.9) to a ridge of high probability regime to one distinct state 
near (—0.45, 0.45), with the change of geometry of potential well. 



and compare them with T 2 from the balanced state S 2 = (—0.45,0.45) (red solid curve) in 



the three-state formulation ( 5TT I. Since the potential well is broader as d increases, the second 



ion (the one which is not trapped in the well) feels its effect more, resulting in an exponential 
increase in the mean escape time. In a channel with descending voltage from left to right 
(U = —0.05 V nm -1 ), it takes a longer time for two ions in the left state (—0.9,0) to escape 
than two ions in the right state (0,0.9), so the black solid curve is above the black dashed 
curve. For two ions in the balanced state — (—0.45,0.45) in the middle of channel, it 
takes an even longer time, so the red solid curve is above the two black curves. 

Notice that the ratio t 2 (-0.9,0)/t 2 (0,0.9) is largest when 0.5 < d < 0.6. Note also that 
for small d the three-state formulation is invalid (so T 2 ( — 0.45,0.45) is meaningless) but that 
for large d the ratio T 2 (— 0.45,0. 45)/t 2 (— 0.9,0) is close to 1: for a broad potential the mean 
escape time is insensitive to the precise initial position of the two ions. 



In Fig. 5.3(b) we plot the left splitting probability p 2 at the left state = (-0.9,0) 



(2) 

(black solid curve) and the right splitting probability 1 — p 2 at the right state 5, = (0,0.9) 
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FIG. 5.3. Fixing U = —0.05 V nm we numerically solve \5.2) and and obtain the mean escape times 

and left splitting probabilities at two states stp = (—0.9,0) and = (0,0.9) for d between 0.1 nm and I nm, 
which determine the escaping rates and transition rates by \i.\2) . We compare it with the results of the three-state 
simplified model |5.1| . (a) The mean escape time ti vs. d at different states, (b) The splitting probability pi vs. d 
at different states, (c) The escaping rates /3' 1 ' = j3 7 ' LI ' and /3' 2 ' = j3, + ' 1,2 '. (d) The transition rates y' 1 ' = y^ 2 ' 1 ' 

and y' 2 ' = yS'' 2 '. Note that those rates only depend on the mean escape time and left splitting probability, and are 
independent of entry rates. 



(black dash-dotted curve) in four-state formulation. For large d, the external potential has 



a large gradient near both ends of the channel (as shown in Fig. 5.1(a) 1, which pulls the 
new ion introduced at either source towards the center of channel, and into balance with the 
existing ion at around (—0.45,0.45). So new entering ions at each source are less likely to 
leave from the same end of the channel, and both the left splitting probability at left state 
(black solid curve) and the right splitting probability at right state (black dash-dotted curve) 
decrease monotonically as d increases. In addition, the left-splitting probability P2 at the 

balanced state si = (—0.45,0.45) in the three-state formulation is plotted by the red solid 
curve, which overlaps with P2(— 0.9,0). As with the escape time, the splitting probability is 
insenstive to the initial position of the ions, and is dominated by the effects of the potential. 

Next we compare the escape rates j3 in the two formulations in Fig. 5.3(c) The black 
solid curve shows the rate at which an ion at left state = ( — 0-9,0) escapes from the 
left side, the black dash-dotted curve plots the rate at which an ion at right state 5, = 
(0,0.9) escapes from right side. The left and right escaping rates at the balanced state 
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,0) 



(—0.45,0.45) in three-state formulation (5.1 1 are plotted by red solid curve and red 



dash-dotted curve, respectively. All escaping rates drop rapidly as d increases. When d > 0.7, 
the potentia l well is so broad that the two ions are trapped in the channel for a long time. In 



i I 1 (2.) 

Fig. 5.3(d) we plot the transition rates between states SV' = (—0.9,0) and 5^ = (0,0.9). 



Due to the inclined voltage, the transition from left state to S^' occurs more often than 



j(2) 



(2 1 i 

the other way, namely the transition rate y\ : ' depicted by solid curve is above J^^' by the 
dash-dotted curve. When d < 0.4, the transition rates are very low (< 10 -4 ns -1 ), so the two 



(1,2) 



states are very distinct. For large d, the observed single state in Fig. 5.2 can be treated as 
average of the two distinct states. 

We remark that the mean escape time of the ion from a one-ion channel Ti (Sj 1 ' = {0}) is 

Q( 1 s ) ti mes larger than T 2 (4' ) X so the esc ape rates fi^ 1 '^ < j3 2 ~ (U) and jS^ 1 ' 1 ' < j3 2 +(U) . 
By ( |3 . 1 3| > we see that for there to be an appreciable probability of having no ion in the channel 



a, 



-a +(u) =A 



- jX. In our example, we choose the smallest 

5(1 



we need %\ (Sj ) _1 

entry rates to be X = jU = 1, so the resulting Pq 1 - is extremely small compared to P^ 1 ' and 

Pq 1 \ and is therefore negligible. Thus, for any entry rates X — 0(1), pi = 0(1), only one-ion 
occupancy (small entry rates) or two-ion occupancy (large entry rates) is observed most of 

time (i.e. P ] 



(i) 



5(l),p(2) 



1). 



In Fig. 5.4(a) fixing X = [X = 5 ns , U = —0.05 V nm 1 , we compare the stationary 
probabilities of two-ion occupancy (solid curves) and one-ion occupancy (dash-dotted curve) 



for various values of d using (i) the four-state formulation (3.8 1 (black curves), (ii) the three- 
state formulation ( |5.1[ ) (red curves), and (iii) by solving the Fokker-Planck equations using 
Cortisol (discrete markers). When d is large, the stationary probabilities obtained from all 
three methods agree with each other, which confirms that the single balanced state in the 
three-state formulation can be treated as an average of the two distinct states (with frequent 
transitions) in the four-state formulation. Because th e mean escape time grows from O(10~ 2 ) 
ns to O(10 ) ns with d increasing as shown in Fig. 5.3(a) a constant entry rate X = jJ. = 5 



ns leads to a transition from initially one-ion dominant to two-ion dominant channel. 

After obtaining the rates and probabilities, we can compare the flux through the channel 
from the rate theory and the solution of the Fokker-Planck equations. We integrate the right 
hand side of ( |2.8b| i with respect to x\ over [—L,L], which yields 

fw-fiL + (n+X)I -(n r^dxi+X [ L Pid Xl )-(f2R-f2L)=0, (5.2) 



-L 



where 



flR 



flL = 




/ dPi ~ d$>i 

D —+Pi - 

\ dxi dx\ 



(—L,x\)dx\ = 



(xi,L)dx\ 



, \ dP\ , , 
(-L)=D-1(-L) 
dx\ 



D^r — (— L,x\)dxi, 
l ax\ 

D^: — (x\,L)dx\. 

L OX 2 



/ dPi ~ d$>\ \ , . dP\ , . 
fiR=D[- r 1 +Pi- r 1 )(L)=D-±(L). 
\ dx\ ax\ I dx\ 



5.1 



In Table 

one by one to the eight terms in (|5.2 



we show that the eight transitions connected to state S^p in Fig. 3.1 correspond 
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Transitions by escaping 


Flux from rate theory 


Flux from Fokker-Planck 


o-(l.l) 

s wh — >s [ l) 

c(l) p l > c(l) 
a l ' ^0 


j3 2 - (M) /f 


/2L 


^+(1.2)^(2) 


-/2« 


P+(l.l)p(l) 


/lL 




-/iff 






M^i 




A/f 


A/j 


-(1,1) 

c(l) "0 „(1) 






+(1,1) 
c(l) a ° , c(l) 




A/ 



Table 5.1 

Transitions computer by the four-state model \3.8\ and the hierarchical Fokker-Planck equation. 




FIG. 5.4. Fixing X = fl =5nj , U = —0.05 V nm~K (a) We plot the s tationary probabilities ofl-ion and 
1-ionby four-state formulation |3.8[ (black curves), three-state formulation ]5.1[ (red curves) and by solving Fokker- 
Planck equation with 28800 elements (discrete marker), (b) We plot the current I vs. d by four-state formulation 
|3.8[ (b lack solid curve) and by three- state formulation \5i\ (red solid curve), which are compared with (//. + fx) 
in \5.3) from solution of the Fokker-Planck equation (discrete marker). 



The left to right flux across the left boundary of the channel is generated by introducing 
new ions /j.Iq + jj, fj^Pi dx\, and the right to left flux across the left boundary of the channel 
is generated by ions leaving the left boundary + fiL- Similarly we have the left to right 
flux across the right boundary of the channel is generated by ions leaving the right boundary 
—/is — f2R, and the right to left flux across the right boundary of the channel is generated by 
introducing new ions A/o + A P\dx\. Thus the overall fluxes are 

ft = iil + At / Pi dxi - fi L - f ZL , f R = -fm - fiR - A/ - A / Pi dxi . (5.3) 



In Fig. 5.4(b) we fix A = \x = 5 ns , U = —0.05 V nm 1 , and plot the current / vs. 
d from the flux (fi +/r) in ( |5.3| l by discrete open diamands, obtained by solving ( |2.8) >. In 
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comparison, the black solid curve depicts the current obtained by applying the four-state rate 
theory ( 3.8 I, and the red solid curve depicts the current by the three-state rate theory ( |5.1| l. We 



see that for a broad potential well with 0.6 < d < 1, all three methods reach a good agreement: 
the four-state works for large d, even though the stationary probability distribution in Fig. 



5.4(a) shows there should be three states for large d. This is because of our use of escape 
times and splitting probabilities to determine the transition rates in the model, rather than by 
trying to estimate hopping rates directly. However, with small d, the three-state model is not 
an accurate description of the Markov process, and the flux quickly becomes inaccurate. 



An important observation from Fig. 5.4(b) is that a maximal current is achieved around 



d = 0.6 nm, which means there exists an optimal shape of the potential well to conduct ions, 
even if the depth of the well remains the same. This result agrees with the argument in 1 1 1 for 
a single ion channel with piecewise linear potential energy. 

We may explain the existence of an op timal flux by looking at the stationary probabilities 
P2 and Pi as a function of d in Fig. 5.4(a) ! When d is small, the potential well is very narrow 



and steep near the binding site, but relatively fiat near the end of channel, so any new ion 
introduced would escape very quickly from the same end by diffusion, leaving the old ion in 
the channel. For example, at d = 0.1 nm, both escaping rates fa {S 2 and fa f^^ 1 ' 2 ^ 

are over 500 ns . Thus the channel has only one ion (Pi > 0.8) most of the time, obviously 
the process of ion entering and leaving from the same end does not generate any through flux. 
On the other hand, when d is large, the potential well is very broad and flat near the binding 
site, two ions can hardly escape from the channel, as shown by the large mean escape time in 



Fig. 5.3(a) once a new ion is introduced to the channel, it quickly moves towards the center, 
and settles into a balanced state in the well with the other ion; thus the 2-ion state dominates 
(P2 > 0.9). In this case the flux is small because two ions are trapped in the channel for a long 
time. 

When neither 2-ion or 1-ion occupancy dominates in the channel, so that there are ade- 
quate transitions between the 2-ion distinct states and frequent escapes from the 2-ion to the 
1-ion state, a large flux is generated. This explains heuristically why an intermediate potential 
well has an optimal geometry. 

Finally we investigate how the entry rates affect the optimal flux for the family of poten- 
tial wells in Fig. J5.1(a)| using the four-state Markov chain formulation ( |3.8| l. Recall that the 
escape rates fa Z^ 1 ' 2 ' and transition rates yP^, ^ are determined by the potential 
through mean escape time and left splitting probability, and thus a re indep endent of the entry 



rates. We fix the applied field U = —0.05 V nm and plot in Fig. 5.5(a) the stationary prob- 
abilities P2 (black) and Pi (red) for A = 1 . 5, 20, 100 ns~ 1 . The intersection points of each pair 
of curves at which P 2 =P\ «0.5 for A = 1,5,20, 100 are at d 0.665,0.6,0.54,0.41 respec- 
tively. This illustrates the fact that when the potential well is narrow and steep at the binding 
site (d small), the escaping rates fa and P 2 are large (shown in Fig. 



5.3(c)!, so the 



entry rates have to increase in order to have equal probabilities of 2-ion and 1-ion occupancy. 
In Fig. 5.5(b) we plot the current 1(d) for A = 1,5,20,100 ns . As expected, the 



current increases as the entry rates increases, but we also find that the value of d at which 
the current is optimised shifts; the critical values of d at which optimal flux is achieved are 
respectively d rj 0.63,0.6,0.59,0.57. Thus the optimal value of d slightly decreases as the 
entry rates increases, which shows that the larger escaping rates of tighter potentials require 
larger entry rates to optimize the flux. 

6. Conclusion. We have presented a set of hierarchical Fokker-Planck equations de- 
scribing ion permeation in multi-ion channels, and reduced these systematically to a discrete 
rate theory. The basis of the reduction is the fact that many channels have internal binding 

23 




FIG. 5.5. Fixing U = —0.05 Vn m . (a) We plot stationary' probability Pi (black) and P\ (red) obtained from 
four-state Markov Chain formulation 3.8i/orA = 1,5,20,100 ns -1 . The intersection points of each pair of curves 
at which P 2 =P 1 ki 0.5 for X = 1,5,20, 100 are at dpi 0.665,0.6,0.54,0.41 nm, respectively, (b) We plot current 
obtained from four-state Markov Chain formulation for p. = 1,5,20,100 ns~ . The critical values of d at which 
optimal flux is achieved are respectively d f» 0.63,0.6,0.59,0.57 nm. 



sites at which ions sit, so that ions transport by hopping between sites on a slow time scale 
while oscillating in the binding sites on a fast time scale. Since the fast oscillation is not key 
in determining the conduction rate, we can reduce the continuous dynamics to the slow tran- 
sition between the discrete states, and thus provide an efficient way to calculate the current 
through the channel. A key component of our reduction was the use of exit times and split- 
ting probabilities to determine the discrete hopping rates, rather than trying to estimate these 
directly using Kramer's theory for example. This means that the predictions of the discrete 
model are accurate even when the internal states are not so well defined. 

In contrast to traditional Eyring rate theory 1111 and the recent study of a one -ion chan- 
nel in HI, we have developed a general theory for multi-ion channels, and have shown an 
intricate coupling between transition rates, mean escape time and splitting probability, due 
to the complexity of the resulting system of Markovian states. The theory is illustrated by a 
two-ion channel, which is the most accessible example that includes the multi-ion complexity. 
We have investigated how conductivity of the channel depends on the diffusion coefficient, 
potential energy landscape, and the ion entry rate. By varying the geometry of the external 
potential while keeping the depth fixed, we observed that when the potential well is narrow 
and steep at the binding site, the 1-ion state dominates, but when it is not the 2-ion state dom- 
inates. In between there is an optimal geometry which maximizes the ion flux by negotiating 
between these two extremes and allowing frequent transitions between the 1-ion and 2-ion 
states. 
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